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ABSTRACT 

We present an extension of the HLLC approximate Riemann solver by Toro, Spruce 
and Speares to the relativistic equations of fluid dynamics. The solver retains the 
simplicity of the original two-wave formulation proposed by Harten, Lax and van Leer 
(HLL) but it restores the missing contact wave in the solution of the Riemann problem. 
The resulting numerical scheme is computationally efficient, robust and positively 
conservative. The performance of the new solver is evaluated through numerical testing 
in one and two dimensions. 
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1 INTRODUCTION 

High energy astrophysical phenomena involve, in many 
cases, relativistic flows, typical examples are superluminal 
motion of relativistic jets in extragalactic radio sources, ac- 
cretion flows around massive compact objects, pulsar winds 
and Gamma Ray Bursts. The modeling of such phenomena 
has prompted the search for efficient and accurate numer- 
ical formulations of the sp ecial relativistic fluid equations 
(for an excellent review see iMartf fc Mul"ieill2003lh There is 
now a strong consensus that the so-called "high-resolution 
shock-capturing" schemes provide the necessary tools in de- 
veloping stable and robust relativistic fluid dynamical codes. 
One of the fundamental ingredients of such schemes is the 
exact or approximate solution to the Riemann problem. 

The solution to the Riemann problem in relativis- 
tic hydrodynamics (RHD henceforth) has been extensively 
studied in literature, and an exact solution can be found 
within high deg r ee of ac curacy by i t erativ e techniques, see 
IMartf fc Miillerl (I1994T) . iPons et alJ (l200(Jt . iRezzolla et"ai] 
(2003) and references therein. One of the major differences 
with the classical counterpart is the velocity coupling intro- 
duced by the Lorentz factor and the coupling of the latter 
with the specific enthalpy. This considerably adds to the 
computational cost, making the use of an exact solver code 
prohibitive in a multidimensional Godunov-type code. 

From this perspective, approximate solvers based on 
alternative strategies have be e n devised: local linea r izatio n 
feulderink fc Mellemal Il995l; iFalle fc Komissarol Il99rj) 
two-s ho ck approximati on jBalsaralll994l :l Dai fc Woodwardl 
Il997l: iM ignone et alJ I2OOFJ) . flux-spli tting methods 
JPonat et alJ Il99gth and so forth, see iMarti fc Miillerl 
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( 2003) for a comprehensive review. Most of these solvers, 
however, rely on rather expensive characteristic decomposi- 
tions of the Jacobian matrix or involve iterative techniques 
to solve highly nonlinear equations. Although they usually 
attain better resolution at discontinuities, some of these 
methods may produce unphysical states wit h negative den- 
sities or pressures, as it has been shown bv lEinfeldt et all 
lll99ll) for linearized Riemann solvers in the context of 
classical hydrodynamics. 

The HLL method devised by iHarten et ail Jl985t for 
classical gasdynamics belongs to a different class of approxi- 
mate Riemann solvers and has gained increasing popularity 
among researchers in the last decade. It has been imple- 
mented in the context of t he relativistic fluid equatio ns by 
ISchneider et~aH <ll993tl and iDuncan fc Hughes! <1994h The 
HLL approach does not require a full characteristic decom- 
position of the equations and is straightforward to imple- 
ment in a Godunov-type code. Besides the computational 
efficiency, this class of solvers has the attractive feature of 
being positively conservative in the sense that preserve ini- 
tially positive densities, energy and pressures. 

The HLL formulation, however, lacks the ability to re- 
solve isolated contact discontinuity and for this reason has a 
more diffusive character than o ther more sophisti cated algo- 
rithms. To compensate for this. lToro et alJ il994T) developed 
an extension of the HLL solver for the Euler equations in- 
troducing a two-state HLL-type solver called HLLC (where 
"C" stands for contact) that improves the treatm ent of the 
contact discontinuity, see also lBatten et all il997f) . Recently 
this approach has been generalized to th e magnetohydrody- 
namic equations <Gurskill20ol IlTI|200FJ| . 

In the present work, we extend this approach to the 
relativistic equation of fluid dynamics. The paper is struc- 
tured as follows: in [J3]the relevant equations are given, in 
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3^1 we describe the new approximate Riemann solver and in 
S^Jwe asses the strength of the new method with one and 
two dimensional tests. 



2 THE RHD EQUATIONS 

The motion of an ideal relativistic fluid is governed by con- 
servation of mass, momentum and energy. The pertaining 
equat ions are cast as a hyper bolic system of conservation 
laws llLandau fc Lifshitdll959l) which, in two dimensions, 
reads 



dU dF x (U) 
dt dx 



dF y (U) 
dy 



0. 



(1) 



where U = (D, m x ,m y , E) is the unknown vector of conser- 
vative variables, whereas F x and F y are, respectively, the 
fluxes along the x and y directions: 



F X (U) 



( Dv x 

rn x v x +j 
m y v x 
, m x 



\ 



) 



( DVy \ 



F y (U) = 



m x v y 

m y Vy + i 



V 



•(2) 



/ 



Generalization to three dimensions is straightforward. 

In equations @, p is the thermal pressure, whereas D, 
m = (m x , m y ) and E are, respectively, the mass, momentum 
and energy densities relative to the lab frame, where the fluid 
has velocity v = (v x ,v y ). Units are conveniently normalized 
so that the speed of light is c = 1. 

The relation between conserved variables U and physi- 
cal quantities V = (p,v x ,v y ,p) is 



D 



IP, 



Dh'jv , E = Dh'y — p , 



(3) 



where p is the proper rest mass density, 7 = (1 — v ■ w)~2 
is the Lorentz factor and h is the specific enthalpy. Proper 
closure is provided by specifying an equation of state in the 
form h = h(p, p). 

For an ideal gas, the enthalpy has the form ph — p + 
f>r/(r — I) and the sound speed is defined by 



(4) 



with r being the (constant) specific heat ratio. By letting 
p/p — > 00, we see that the square of the sound spee d has 
the limiting value cj — > T — 1. Since it can be shown jTaubl 
ll948l:lAmle ll989HMignone et all20oH) that the specific heat 
ratio r cannot exceed 2, one always has (? s < 1. This is an 
important result for the positivity of the HLL and HLLC 
schemes and will be used in a later section. 

Equation gives U in terms of the primitive state 
vector V . The inverse relation involves the solution of a 
nonlinear equation for the pressure p: 



E + p = Z>y 



rP7 



(•») 



where 7 = [l — \m\ 2 /(E + p) 2 ] 2 ■ Equation can be 
solved by any standard root finding algorithm. 




Figure 1. Graphical representation of the Riemann fan in the 
x — t plane. The two initial states Ul and Un decay into two 
nonlinear waves (with speeds \l and and a linear contact 
wave with velocity A*. The resulting wave pattern divides the 
x — t plane into four regions each defining a constant state: Ul , 
U* L , U* R and Ur. 



2.1 The Riemann Problem in RHD 

Consider a conservative discretization of £Q| along the x- 
direction: 



■W 



fi 



f, 



(6) 



At n Axi 
The numerical flux functions f i+ i follow from the so- 
lution of Riemann problems with initial data: 



U 



U(x,Q) 



U 



R,i+i 



if x < x, 
if x > X- 



(7) 



and U R i+ i are the left and right edge values 



where U L 
at zone interfaces. 

The solution of the Riemann problem for the spe- 
cial relativistic fluid eq uations has been i nvestigated by 
iMarti fc MuUeil dl994l) . iPons et alJ feOOOh . iRezzolla et ail 
It consists of a self-similar three- wave pattern gener- 
ated by the decay of the initial discontinuity 0. The result- 
ing Riemann fan (Fig.0 is bounded by two nonlinear waves 
(representing either shocks or rarefactions) separated by a 
contact discontinuity moving at the fluid velocity. Across 
the contact discontinuity, pressure and normal velocity are 
continuous whereas density and tangential velocities expe- 
rience jumps. The same holds also in the non-relativistic 
limit. Contrary to the Newtonian counterpart, however, all 
variables are discontinuous across a shock wave or c hange 
smoothly through a rarefaction fan JPons et alJl200d) . This 
is a consequence of the velocity coupling introduced by the 
Lorentz factor 7 and by the coupling of the latter with the 
specific enthalpy h. 

The resulting wave pattern can be solved to a high de- 
gree of precision by iterative techniques and has been imple- 
mented for th e first time in th e one-d imensional Godunov- 
type code by IMarti' fc Miillerl dl99rl) . Nevertheless, when 
tangential velocities are included, the computational effort 
increases considerably, and the use of an exact solver in 
a multidimensional Godunov-type code can become pro- 
hibitive. 

Here we consider a diff erent approach, bas ed on the 
original prescription given bv lHarten etafl Jl983l) (HLL) for 
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the classical Eule r equations and subsequently extended by 
iToro et alJ il994h (HLLC). Both the HLL and HLLC formu- 
lations do not require a field by field decomposition of the 
relativistic equations, a feature which makes them particu- 
larly attractive, specially in multi-dimensional applications. 



3.1 HLLC Solver 

The HLLC scheme restores the full wave structure inside the 
Riemann fan by replacing the single averaged state defined 
by © with two approximate states, U* L and U* R . These 
two states are separated by a middle contact wave which is 
assumed to have constant speed A* , so that the full solution 
to the Riemann problem now reads 



3 THE HLL FRAMEWORK 



Harten, Lax and van Leer (lHarten et alJll983l) proposed an 
approximate solution to the Riemann problem where the 
two states bounded by the two acoustic waves are averaged 
into a single constant state. In other words, the solution to 
the Riemann problem on the x/t = axis consists of the 
three possible constant states: 



U(0,t) 



where we dropped, for sim plicity, the half integer notation 
i + jr. lHarten et alJ <1983l) noted that the single state U h " 
could be constructed from an a priori estimate of the fastest 
and slowest signal velocities Xl and X R : 



C7(0, t) 



Ul 

u R 

Ur 



if 
if 
if 
if 



Xt > o, 

X L < s= A* , 
A* < s= X R , 
Xr < 0, 



(12) 



and the corresponding intercell numerical fluxes are: 



if 



X L > 0, 





if 


X L > 0, 






n 


if Xl < < A* 


jj-hll 


if 


Xl < < Xr , 


(8) 


/= < 


F" R 


if A* < < Xr 


Ur 


if 


Xr : SC , 






Fr 


if Xr < . 



(13) 



The intermediate state fluxes F* L and F* R may be ex- 
pressed in terms of U* L and U* R from the Rankine-Hugoniot 
jump conditions: 



A (17* — U) = F* — F , 



(14) 



jjhll _ XrUr 



XlUl+F l 



Xr — Xl 



(9) 



where F L = F X (U L ), F R = F X (U R ). Notice that equation 
© represents the integral average of the solu tion of the 
Riemann problem over the wave fan (lTorolll997l) . 

The corresponding interface numerical flux is defined 

as: 



where 



F L 


if 


X l ^0, 




jphll 


if 


Xl^O^Xr, 


(10) 


Fr 


if 


Xr^O, 




XrFl 


-Xl 


Fr + XrX l {Ur~Ul) 


(11) 






Xr — Xl 



where here and throughout the following, quantities without 
a suffix L or R refer indifferently to the left (L) or right 
(R) states. Notice that, in general, F = F X (U) but F* / 
F X (U*). 

We remind the reader that the HLL and HLLC solvers 
differ in the representation of the intermediate states. In the 
case of supersonic flows (Xl > or A_r < 0), in fact, the two 
solvers become equivalent. The same result also holds for an 
exact Riemann solver. 

If Xl and Xr are given (see i!3.1.1l . equation 1141 rep- 
resent a system of 2n equations (where n is the number of 
components of U) for the 4n+l unknowns U* L , U* R , F* L , F* R 
and A* . Three additional constraints come from the require- 
ments that both pressure and normal velocity be continuous 



across the contact wave (i.e. 



• P*r — Pl) an( i that 



Thus, given a wave speed estimate for the fastest and 
slowest speeds Xr and Xl (see H3.1,ip . an approximate so- 
lution to the Riemann problem can be constructed and the 
intercell numerical fluxes for the conservative update © are 
computed using (1101 . 

This approach has been applied for the first time to the 
one-di mensional re lativistic equations bv | Schneider et alJ 
( 1993) and later bv lDuncan fc Huehesl dl994T) for the multi- 
dimensional case. 

Although the HLL prescription is computationally in- 
expensive and straightforward to implement, a major draw- 
back is its inability to resolve contact or tangential waves. 
On the con t rary t he HLLC scheme, originally introduced by 
IToro et alJ dl994h in the context of the Euler equations of 
classical gasdynamics, does not suffer from this loss. In the 
next section we generalize this approach to the equations of 
relativistic hydrodynamics. 



A* = v* L = Vx,R- This, however, yields a total of only 2n + 3 
equations, still not sufficient to solve the system. In order 
to reduce the number of unknowns and have a well-posed 
problem, further assumptions have to be made on the form 
of the fluxes F* . Here we assume that the two-dimensional 
fluxes can be written as 



\ 



D*v* 



(15) 



In such a way, both U* and F* are expressed in terms of 
the five unknowns D* , v*, m*, E* and p*. The normal com- 
ponents of momentum in the star region, m* L an d rn^, R , are 
not independent variables since, for consistency, we require 
that m* = {E* +p*)Vx- In the classical case, this assumption 
becomes equivalent to m* = p*X* . This yields a total of 11 
equations in 11 unknowns. 
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Writing explicitly equation 1141 for the left or the right 
state yields 

D*{A~A") = D(A — v x ) , 

m*(A - A*) = m x (A — v x ) + p* — p, 

mJ(A-A*) = my(\ — v x ), 

£T(A-A*) = E(\-v x )+p*\* -pv x . 

If one combines the last of 1 Kit together with the second 
one, the following expression giving A* in terms of p* may 
be obtained: 



(16) 



(A-Ap*)v* = B+p* 



(17) 



where A = AE — m x , B — m x (A — v x ) — p. 

By imposing p* L = p* R across the contact discontinu- 
ity we find the following quadratic equation for A*: 



nhll / \ * \ 2 / -w—^hll , TnhXX \ \ * i fill r\ 

F E (A) -{E + F m jA + m x =0. 



(18) 



In equation 1181 . and F„J are the energy and mo- 
mentum components of the HLL flux given by equation 1111 . 
whereas E hu and m x 11 are the energy and normal momen- 
tum components of the HLL state vector, equation Of 
the two roots of equation 1181 only the one with the mi- 
nus sign is physically acceptable, since it lies in the range 
(— 1, 1) and, according to the wave-speed estimate presented 
in can be interpreted as an average velocity between 

Al and Ar. The mathematical proof of this statement is 
given in the appendix 0. Besides, the same root has the 
correct classical limit, that is A* — » m£ /p as v/c — > 
and h — > 1 . This wave speed is the same one proposed by 
iTorol (Il997f) and further discussed in lBatten et al] (Il997l) . 

Once A* is known, p* is computed from 1171 and the 
components of U* are readily obtained from 1161 . 

Finally we notice that the method is consistent, in that 
the integral average over the Riemann fan auto matically sat - 
isfies the consistency condition by construction ||ToroJl997): 



(a* — \l)u* l + {\r — a*)u* r = uhu 

Ar — Al 
or, alternatively, 

F* l Ar(A* - A L ) + F R A L (A R - A* 



Ar — At 



= A*F 



(19) 



(20) 



Incidentally, we notice that equation I18H could have been 
obtained by algebraic manipulations of equations 119H and 



3.1.1 Wave Speed Estimate 

The wave speeds needed in our formulation are estimates 
for the lower and upper bounds of t he signal ve locities in 
the solution to the Riemann problem ]Torolll9'9^1 . Here we 
co nsider the rel ativistic generalization of the estimates given 
bv lDavis! 1^988) for the Euler equation of g asdynamics. The 
same c hoice has been initially ad opted bv ISchneider et all 
lll993h . iDuncan fc Hughes! ]1994T) in their relativistic HLL 
solver and is commonly used by other authors, see, for ex- 
ample, [De^anm^^^uccimitiifil l|2002|l. Specifically we set: 

Al =min(A-(V H ),A-(Vi)), 

(21) 

A R = max(A + (V fl ),A + (Vi)), 



where A+ and A_ are the maximum and minimum eigenval- 
ues of the Jacobian matrix dF/dU. They are the roots of 
the quadratic equation 



(A — v x ) 2 = <j s (l - A 2 ) , 

with a s = c 2 /(7 2 (l — c 2 )), and hence 



A±(V) 



1 +cr s 



(22) 



(23) 



It should be mentioned that the wave speed estimate 
1211 is not the only possible one and different choices (such 
as the Roe average) may be considered. 



3.1.2 Positivity of the HLLC scheme 

Adopting the same notations as in lBatten et al] dl997h . we 
denote with G the set of physically admissible conservative 
states: 



{ ( D \ 



G= { 



D>0 



E > y/m x + m 2 y + D 2 



(24) 



where the second inequality simultaneously guarantees pres- 
sure positivity and that the total velocity never exceeds the 
speed of light. 

We remind the reader that the pressure p(U*) com- 
puted from the conservative state U* using Jjy should not 
be confused with p* appearing in the flux definition 115H . 
The two pressures are, in fact, different and the positiv- 
ity argument should apply to p(U*) rather than p* , which 
can take negative values under certain circumstances. Sim- 
ilar considerations hold for the velocity A* of the contact 
discontinuity for which, in general, we have A* 7^ v x (U*). 
Thus p* and A* may be more conveniently considered as 
auxiliary variables. 

This is one of the fundamental differences between our 
relativistic solver and the classical HLLC scheme, for which 
A* = m x I p and thus only p* plays the role of an auxil- 
iary parameter. This behavior is a direct consequence of the 
relativistic coupling between thermodynamical and kinetical 
terms, a feature absent in the Newtonian formulation. 

The positivity of the HLLC scheme is preserved if each 
of the two intermediate states U* L and U* R are contained in 
G. 

For the density, the proof is trivial and follows from the 
inequalities Al ^ A* ^ Ar and Al ^ v x (L,R) ^ Ar, see 
Appendix 3Al 

Unfortunately the analytical proof of the second state- 
ment presents some algebraic difficulties, since the second 
of 12411 reduces to an inequality for a quartic equation in 
A*. However, extensive numerical testing, part of which is 
presented in has shown that the second of 1241 is al- 
ways satisfied for all pair of states Ul and Ur whose wave 
speeds are computed according to 1211 and for which an ex- 
act analytical solution to the Riemann problem exists (i.e. 
no vacuum is created). 



An HLLC Riemann Solver for Relativistic Flows: I. Hydrodynamics 5 



4 ALGORITHM VALIDATION 

We now provide some numerical examples to test our new 
HLLC solver. For th e test problems considered i n this sec- 
tion we closely follow ILucas-Serrano et all i2004f) . 



4.1 Implementation Details 

The numerical integration of the relativistic equations £Q| 
proceeds via the conservative update (|HJ . For the first-order 
HLLC scheme, we compute the inter-cell numerical fluxes 
using l|10|l with left and right states given, respectively, 
by Ui and C/i+i- 

For the second order scheme, the input to the Riemann 
problem are the states 

8V7 



y n +? _ y' 



: _ y" r 2 
,R *+l 



+ 



2 ' 



(25) 



where V i 2 follows from a simple Hancock predictor step, 
At" 



2Axi 



F V 



— F [V 



(26) 



with V n i and V™ 1 computed from 1251 1 by replacing 

*~ p ^ , L; t — ^ , IX 

V n+ i with V n . 

The SV's appearing in equation 1251 are computed at 
the be ginning of the time step using t he fourth-order limited 
slopes <IColellalll985l ISaltzmarJI 19941 : 



SVi 



-A V t 



SVi+i + SVi-t 



, AiVi 



where 

A,V i = amin(|AV 4 |,|AV i _i|) , 
and SVi are the second-order slopes 
SVi = s I min(A i y i ,jA y s |) , 

AF» = V %+1 - V, , A Vi 



Vi+i - V t - 



sign(AVi) + sign(AVi-i) 



(27) 

(28) 

(29) 
(30) 

(31) 



The parameter a £ [1,2] adjusts the limiter compression, 
with a = 2 (a = 1) yielding a more (less) compressive lim- 
iter. Notice that, although the use of fourth-order slopes at- 
tains sharper representations of discontinuities, the scheme 
retains global second-order spatial accuracy. 

We do not make use of any artificial steepen- 
ing al gorithm to enhance re solution across a contact 
wave llMarti' fc Mulled Il996t ILucas-Serrano et al.1 12004 * 



iMignone et alJl2005T) in order to highlight the intrinsic ca- 
pabilities of our new HLLC solver. In the one-dimensional 
tests, the computational domain is the interval [0, 1] and the 
compression parameter is a = 2. In two dimensions we set 
a = 2, q = 1.25 and a = 1 for density, velocities and pres- 
sur e, respectively. A dditional shock flattening, computed as 
in jMartf fc Mulledll99rJ) . is used in Sg^land 34lflto pre- 
vent spurious numerical oscillations. Outflow boundary con- 
ditions are set in problem 1-4. 
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Figure 2. Comparison between the HLL (dotted line) and HLLC 
(dashed line) Riemann solvers for problem 1 at t = 0.4. Only the 
density profiles are shown. Computations were performed with 
the first-order scheme on 100 computational zones with CFL = 
0.8. The solid line gives the analytical solution. The major differ- 
ence between the two approaches is the resolution of the contact 
wave. 



Multidimensio nal integratio n is achieved via Strang di- 
rectional splitting jStrandll968lh that is, by successively ap- 
plying one-dimensional operators in reverse order from one 
time step to the next one, i.e. U n+2 = C x C v U n+1 and 
U n+1 — Ly£. x U n . Here £. x is the operator corresponding 
to the conservative update © (and similarly for C y ). The 
same time increment At should be used for two consecutive 
time steps. 



4.2 Problem 1 

The first test consists in a Riemann problem with initial 
data 



(p,v x ,p) = 



(1,0.9, l) for x<0.5, 
(1,0, 10) for x>0.5. 



(32) 



Integration is carried with CFL = 0.8 until t — 0.4 and an 
ideal equation of state with T = 4/3 is used. The breakup of 
the discontinuity results in the formation of two shock waves 
separated by a contact discontinuity. 

In Fig. |21 we plot the analytical solution for the rest 
mass density together with the profiles obtained with the 
first-order HLLC and HLL schemes on 100 uniform compu- 
tational zones. The two integrations behave similarly near 
the shock waves, but differ in the ability to resolve the con- 
tact discontinuity. As expected, the HLLC scheme yields a 
sharper representation of the latter, whereas the HLL solver 
retains a more diffusive character. 

The Li norm errors of density are shown in the top- 
left panel of Fig. |H] for different resolutions. For the sake of 
comparison, computations have also been performed using 
the more sophisticated exact Riemann solver describ ed in 
the one-dimensional code bv iMartf fc Mulled (Il996lh The 
errors obtained with the present HLLC scheme and the exact 
Riemann solver are comparable at low resolution (~ 15.3% 
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0.0 0.2 0.4 0.6 



Figure 3. Numerical solutions obtained with the second-order 
HLLC scheme applied to problem 1. The solid line represents the 
analytical solution, while computed profiles of density (crosses), 
pressure (plus signs) and velocity (filled circles) are shown on 400 
computational zones at t = 0.4. A CFL number of 0.8 was used. 



Figure 4. Computed density profiles for the first-order HLL 
(dotted line) and HLLC (dashed line) schemes for problem 2 at 
t = 0.4. We used 100 grid points and CFL number of 0.8. The an- 
alytical solution is plotted as a solid line. As expected, the HLLC 
solver behaves quantitatively better than the HLL scheme across 
the contact wave. 



and ~ 13.6% respectively on 100 points) and become nearly 
identical as the number of points increases. Conversely, the 
errors computed with the relativistic HLL scheme are bigger 
(~ 22.2% on 100 points) and show that almost twice the 
resolution is needed to achieve the same accuracy obtained 
with the HLLC or the exact solver. 

Fig.[3]shows the results computed with the second-order 
HLLC scheme on 400 zones, at the same time. The exact 
profiles for density, velocity and p ressure are plotted as solid 
lines. Additional slope flattening (iMartf fc M iillcr 1996J has 
been used to reduce the spurious numerical oscillations ob- 
served behind the shock front. All discontinuities are ade- 
quately captured and resolved on few computational cells, 
~ 3 for the shocks and ~ 4 — 5 for the contact discontinuity 
(contrary to ~ 7 when the HLL solver is employed). 

The error in Li norm is ~ 2.3% for 400 grid zones and it 
has been computed at different resolutions using the HLL, 
HLLC and exact Riemann solver, see Fig. |5] Not surpris- 
ingly, the second-order interpolation considerably reduces 
the errors and higher convergence rates are expected for all 
schemes. Nevertheless, the three solvers mostly differ in the 
resolution at the contact discontinuity and, for n ^ 800 grid 
points, the HLLC and exact Riemann are practically in- 
distinguishable, while at the maximum resolution employed 
(3200 zones), the error computed with the HLL scheme is 
still ~ 20% bigger. 



4.3 Problem 2 

In the second test, we prescribe the initial condition 
(l, -0.6, 10) for x < 0.5 , 
( 10, 0.5, 20) for x > 0.5 , 



(33) 




with an ideal equation of state with F — 5/3. Integration 
stops at t — 0.4 and CFL = 0.8 has been used in the inte- 



0.0 0.2 0.4 0.6 0.8 1.0 



Figure 5. Second-order HLLC scheme applied to problem 2 at 
t = 0.4 on 400 computational zones with CFL = 0.8. As in Fig. [3] 
profiles for density, pressure and velocity are plotted using crosses, 
plus signs and filled circles. The contact wave is the only discon- 
tinuity in the solution and is clearly visible at x m 0.4. 



gration. The initial discontinuity evolves into left-going and 
right-going rarefaction waves with a contact discontinuity in 
the middle. 

Results for the first-order HLL and HLLC schemes on 
a 100-point uniform grid are shown in Fig. [I] Again, notice 
the sharper resolution of the HLLC scheme in proximity of 
the contact wave. The smooth rarefaction waves are equally 
resolved by both schemes. 

The behavior of the solution under grid resolution ef- 
fects is described in the top-right panel of Fig. |H| Since the 
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Figure 6. Results of the second-order HLLC scheme applied to 
shock tube problem 3 at t = 0.4 on 400 computational zones. 
Integration has been carried with CFL = 0.8. The solution is 
comprised of a left-going rarefaction wave, a right-going contact 
and shock wave moving close to each other. 



- p/1000 


- 






/ X > 






< 


-P/10 / >w 

/^^^^— i^^^ M ____^****Htt»«»...^ - > 




1 i i r I r i i 1 i i t 1 i 





0.0 0.2 0.4 0.6 0.8 1.0 



x 

Figure 7. Computed profiles of density, pressure and velocity for 
the second-order HLLC scheme applied to blast wave problem 
(problem 4). Integration has been carried with CFL = 0.8 on 
400 uniform zones until t = 0.4. The configuration is similar to 
that of problem 3, but the shock and the contact waves are now 
much closer to each other. 



only discontinuity in the problem is the contact wave, the L\ 
norm reflects mostly the different resolution across the dis- 
continuity. The HLLC and the exact solver perform nearly 
identically, while the HLL exhibits a slightly slower conver- 
gence rate. At the maximum resolution, the error in the HLL 
scheme is ~ 4.3% to be compared to the ~ 3.0% and 3.1% 
errors obtained from the other two Riemann solvers. 

These differences are again reduced in the second-order 
HLLC scheme, Fig. |SJ for which the convergence rates are 
similar, as shown in the top-right panel of Fig. [5] 



find strong differences between our HLLC method and the 
HLL scheme. Resolution effects are given in the bottom left 
panels of Fig. |H] and |§] for the first-order and second order 
schemes, respectively. As one can see, the solutions com- 
puted with the HLL, HLLC and exact solvers behave nearly 
in the same way, with the L\ norm errors being different by 
less than 1% at low resolution and becoming identical for 
n ^ 800 grid points. We believe that this might be due to 
the proximity of the contact and shock waves. The quality 
of our results is, however, similar and comparable to those 
obtained in previous studies. 



4.4 Problem 3 



4.5 Problem 4 



The initial condition for this test is 

(10,0,40/3) for x<0.5, 
(1,0,0) for x>0.5, 



(34) 



with r = 5/3. For numerical reasons, the pressure in the 
left state has been set equal to a small value, p — 2/3 • 10 -6 . 
Integration is carried with CFL = 0.8 on 400 grid points; 
the final integration time is t = 0.4. The initial configuration 
results in a mildly relativistic blast wave, with a maximum 
Lorentz factor 7 max ~ 1.4. The Riemann fan consists of a 
rarefaction wave moving to the left, a shock wave adjacent 
to a contact discontinuity, both moving to the right, see Fig. 

m 

Our HLLC scheme is able to capture discontinuities 
properly; in particular, the shock is resolved within 2-3 zones 
and the contact discontinuity smears out over 4-5 zones. We 
remind again that the interpolation algorithm does not make 
use of additional artificial compres sion to enhance resolu- 
tion across the contact wave, as in iMartf fc Mulled <1996f) . 
Moreover, we repeated the test also with the exact Riemann 
solver and did not find any noticeable difference. In addition 
and contrary to the previous two test problems, we did not 



(p,v x ,p) 



(35) 



In the fourth shock-tube, we prescribe the following initial 
discontinuity 

(1,0, 10 3 ) for x<0.5, 
(1,0, 10 -2 ) for x>0.5. 

Again we adopt an ideal equation of state with V = 5/3. 
The resulting pattern is similar to that of problem 3, but 
the specific enthalpy in the left state is much greater than 
unity, thus resulting in a more thermodynamically relativis- 
tic configuration. The solution computed with the second- 
order scheme at t = 0.4 is shown in Fig. Qon 400 computa- 
tional zones and CFL=0.8. 

The high pressure jump produces a strong shock wave 
and a contact discontinuity very close to each other, moving 
to the right at almost the same speeds. The higher com- 
pression in the shell is due to a more pronounced relativistic 
length-contraction effect caused by a higher Lorentz factor, 
7max ~ 3.7. The smaller thickness of the shell between the 
shock and the contact wave makes this test numerically chal- 
lenging and particularly demanding in terms of resolution. 

Our relativistic HLLC scheme is able to reproduce the 
solution within a satisfactory agreement, even without us- 
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Figure 8. Grid resolution studies for the first-order schemes us- 
ing the HLL (filled circles), HLLC (crosses) and exact Riemann 
solvers (triangles). Results are shown for problem 1 (PI, top left), 
problem 2 (P2 top right), problem 3 (P3, bottom left) and prob- 
lem 4 (P4, bottom right). Computations have been performed on 
50, 100, 200, 400, 800, 1600 and 3200 grid zones with the same 
CFL number (0.8) for all runs. Notice that both the HLLC and 
exact solvers perform better than the HLL scheme in problem 1 
and 2, while all schemes yield nearly identical results in problem 
3 and 4. 



ing a contact steepening algorithm. The absolute global er- 
ror in density is 6.5% and the density peak in the thin shell 
achieves ~ 81.6% of the exact value. Our results are there- 
fore similar to previous ones proposed in literature. 

It should also be pointed out that, similarly to prob- 
lem 3, we did not find any improvement in the solution by 
switching to the exact Riemann solver or using the relativis- 
ts HLL scheme. This is confirmed by the resolution study 
carried for the first and second-order schemes (bottom right 
panels in Fig.lHJandEJ. Again, we suggest that the ability to 
capture the discontinuities relies on the interpolation prop- 
erties of the algorithm and has a weaker dependence on the 
Riemann solver for this particular class of problems. 

4.6 Relativistic Planar Shock Reflection 

The initial configuration for this test problem consists in a 
cold (p — 0), uniform (p = 1) flow impinging on a wall lo- 
cated at x — 0. The flow has constant inflow velocity Vi n and 
the reflection results in the formation of a strong shock wave. 
For t > the shock propagate s upstream and the s olution 
has an analytic form given by iMartf fc Mull erl ll996h : 



p(r,t) = 



for r > v s t . 
for r < v a t . 



where 

r + i + r ( 7in - i) 

a = 

r - 1 



(r-i 



'yin fain \ 



(36) 



(37) 



7in + 1 

are the compression ratio and the shock velocity, respec- 
tively. Behind the shock wave (r < v s t), the gas is at 



Figure 9. Grid resolution studies for the second-order schemes. 
The notations are the same ones used in Fig. [H] Computations 
were obtained using the same CFL number (0.8) for all cases. 




Figure 10. Relativistic planar shock reflection test with CFL = 
0.4 and 100 grid points. Results are shown at t = 1.5 for the 
second-order HLLC scheme. The cold supersonic flow enters at 
x = 1 and a reflective boundary condition is imposed on the left, 
at x = 0. The reflected shock is located at x = 0.5. 



rest (i.e. v — 0) and the pressure has the constant value 
p(r,t)(-fin — l)(r — 1). Conversely, in front of the shock all 
of the energy is kinetic and thus p = 0, V — Vin . 

For numerical reasons, pressure has been initialized to a 
small finite value, p = e(T- 1), with e = 10~ 10 and F = 4/3. 
The computational domain is covered with 100 computa- 
tional zones and the initial inflow velocity is Vi n = —0.99999 
corresponding to a Lorentz factor fi„ ~ 224. Integration is 
carried with CFL =0.4. 

Fig. 1101 shows the solution at t — 1.5, after the shock 
has propagated Ax s « 0.5 from the wall. The relative 
global errors (defined as e(Li)/ £\ p ex (xi)Axi) for den- 
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sity, velocity and pressure are, respectively, 1.8%, 1.4% 
and 1.4%. This result is in excellent quantitative agree- 
ment with the numerical so l utions obtained by other au- 
thors jMarauina et aT]|l992t iMartf et alJll997l : lAlov et all 
ll999l:lDel Zanna fc Bucciantinill2002l) . In this test, similarly 
to problem 1, shock flattening was employed to prevent nu- 
merical oscillations. 

Also, from the same figure, we notice that our solver suf- 
fers from the wall-heating phenomenon, a common pathol- 
ogy in modern Godunov-type schemes. The degree of this 
pathology is higher than the HLL scheme but less than the 
exact Riemann solver. We also point out that the problem 
may be partially mitigated by a proper fine-tunings of the 
parameters involved in the reconstruction and steepening al- 
gorithms. However, we did not follow that approach in the 
present work. 



4.7 Two-Dimensional Riemann Problem 

Two dimensional Riemann problems involve the interactions 
of four elementary waves (either shocks, rarefactions, and 
contact discontinuities) initially separating four constant 



states. They hav e been formulated bv ISchulz et alJ lll993ft . 
lLax fc Liul (Il998l) in the context of classical hydrodynamics. 
He re we consider a relativist i c exte nsion, initially proposed 
bv lDel Zanna fc Bucciantinil (120021) . where the initial con- 
figuration involves two shocks and two tangential disconti- 
nuities. 

The domain is the square [—1, 1] x [—1, 1] covered with 
400 2 computational zones. The four quadrants NE (x, y > 
0), NW (x<0< y), SW (x,y< 0), SE (y < < x) divide 
the square into four constant-state regions: 



(p,v x ,v v ,p) 



(0.1,0,0,0.01) for x,y>0, 

(0.1,0.99,0,1) for x<0<y, 

(0.5,0,0,1) for x,y<0, 

(0.1,0,0.99,1) for y<0<x. 



(38) 



We use an ideal equation of state with T = 5/3. The inte- 
gration is carried out with CFL=0.4 till t = 0.8. 

Notice that the initial condition 1381 1 does not exactly 
prescribe two simple shock waves at the NW — NE and 
SE — NE interface. The correct version of this problem has 
been considered by iMignone et alJ ll2005l) . For the sake of 
comparison, h owever, we chose to adop t the same initial 
condition as in lDel Zanna fc Bucciantinil J2002TI . 

The top and bottom panels in Fig. llll show. respectively, 
the solutions computed with the HLLC and HLL solvers. 
The breakup of the initial discontinuity results in two equal- 
strength curved shock fronts propagating from regions NW 
and SE into the upper right portion of the domain (NE), 
top panel of Fig. ^3 Region SW is bounded by two tan- 
gential discontinuities and a jet-like structure emerges along 
the main diagonal. 

The initial steady tangential discontinuities, located at 
the W and S interfaces, remain automatically sharp in the 
HLLC formulation, since they are exactly captured by the 
approxim ate Riemann so l ver. T he same results has also been 
shown bv lMignone et al.l (120051) who used a two-shock itera- 
tive nonlinear solver. We emphasize that this is property per- 



0.5 



0.0 





Figure 11. Density logarithms for the two-dimensional Riemann 
problem on 400 2 zones at t = 0.8; the top (bottom) panel shows 
the results obtained with the second-order HLLC (HLL) scheme. 
Thirty equally spaced contours are shown. Integration has been 
performed with CFL = 0.4. Curved transmitted shocks are vis- 
ible in the upper right portion of the domain. The drop-shaped 
region in the lower left portion is bounded by two tangential dis- 
continuities. The HLL (bottom panel) shows additional numerical 
diffusion which degenerate into spurious waves propagating along 
the main axis. 



tains to the Riemann solver itself and does not depend on the 
interpolation algorithm. Indeed the same result holds when 
the first-order scheme is employed. This feature is absent 
from the HLL formulation, where tangential discontinuities 
spread along the cartesian axis due to extra numerical diffu- 
sion. This behavior is manifestly evident in the grid-aligned 
spur ious waves visible in the b ottom panel of F ig. fTTI see 
also | Del Zanna fc Bucciantinil (120021) ; iLucas-Serrano et alJ 
(2004). 
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4.8 Axisymmetric Relativistic Jet 

Finally, as an example of an astrophysical application, 
we consider the propagation of a light, axisymmetric rel- 
ativistic jet in 2-D cylindrical geometry. For the sake 
of comparison, the parameters of the simulatio n are the 
same used by |Del_Za nna fc Bucciantinil i2002f) and by 
iLucas-Serrano e"t^dT ~ il2004i) . The initial condition is pre- 
scribed as 

{(O.l, 0, 0.99, 10~ 2 ) for r,z <1, 
(10, 0, 0, 10- 2 ) otherwise. (39) 

The jet is pressure matched and its internal relativistic Mach 
number is 17.1. We use an ideal equation of state with 
r = 5/3 both for the jet and the ambient medium. The com- 
putational domain covers the region $J r 12, ^ z ^ 35, 
with 240 x 700 grid points, so that we have 20 cells per jet 
radius. At the symmetry axis, r — 0, we impose reflecting 
boundary conditions; outflow boundary conditions are set 
everywhere else, except at the inlet region where we keep 
the beam constant. The CFL number is 0.5 and the jet evo- 
lution is followed until t = 80. 

For the sake of comparison, we have also carried the 
simulation using the relativistic HLL solver. Figure lT2l shows 
two snapshots of the rest mass density at times t — 40 and 
t — 80. The upper-half of each panel refers to the HLLC 
integration, whereas the lower portion displays the result 
obtained with the HLL scheme. We see that all the struc- 
tural features characteristic of jet propagation can be clearly 
identified, with good resolution of shock waves and contact 
discontinuities. It is clear from the figures that the HLLC in- 
tegration shows a significantly greater amount of small scale 
structure, that is not visible in the HLL results. This is due 
to the larger numerical diffusion introduced by the latter in 
subsonic regions which prevent sharp resolution of shear and 
tangential waves. 

The average advance speed of the jet head is ~ 0.39 
(to be c ompared with a on edimensional theoretical estimate 
of 0.44: iMartf et alJ Jl997ft '). Moreover, we can observe the 
absence of the carbuncle problem, that usually appears as a n 
extended nose in front of the jet, on the axis iQuir ki ll 9941. 



5 EFFICIENCY COMPARISON 

Previous numerical tests have shown that the quality of 
solution achieved with the HLLC scheme can be com- 
petitive with more complex exact or iterative non-linea r 
Riemann solvers, see for example IMartf fc Mulled |l996), 
iMignone et alJ <l2005ll . Another aspect which plays in favour 
of the HLLC formalism is the computational efficiency, par- 
ticularly crucial in long-term simulations in two and three 
dimensions. 

Table^gives the normalized CPU time required by the 
HLL, HLLC and approxim ate two-shock nonline ar Riemann 
solvers (for the latter see IMignone et alJ [2005) . All three 
solvers are available in the C author's code and have been 
written performing similar degree of optimizations. On the 
contrary, the FORTRAN code for the exact solutio n to th e 
Riemann problem, available from IMartf fc Mulled {2003), 
was found to be more than a factor of 7 slower than the HLL 
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Figure 12. Gray scale images of the rest-mass density logarithm 
for the relativistic jet simulation at t = 40 (top panel) and t = 80 
(bottom panel). In each panel, the HLLC (HLL) solver has been 
used for the upper (lower) portion of the image. The resolution is 
20 points per jet radius, corresponding to grid size of 240 X 700 
computational zones. Integration has been carried with CFL = 
0.5 



solver. We believe that this might be due to a lower degree of 
optimization and by the time c onsuming numeri cal integra- 
tion across the rarefaction fan iPons et alj EoOOI. For illus- 
trative purposes, we consider the first four one-dimensional 
tests and the two-dimensional Riemann problem. Integra- 
tions have been carried using the first-order scheme on 4000 
and 400 2 zones, respectively. No optimization flags were used 
during the compilation. 

From the table, one can easily conclude that the HLLC 
scheme requires little additional costs with respect to the 
HLL approach (between 4% and 7% in one dimension) , while 
the iterative nonlinear solver is certainly more expensive, 
being by a factor of more than 30% slower. 

In making the comparison, however, one should keep 
in mind that HLL-type solvers are iteration-free since the 
underlying algorithms always require a fixed number of op- 
erations, regardless of the initial condition. On the contrary, 
iterative nonlinear Riemann solvers have a certain degree of 
adaptability since the number of iterations to achieve con- 
vergence depends on the strength of the discontinuity at the 
zone interface. In smooth regions of the flow, for example, 
fewer iterations are usually needed. This explains why the 
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Test 


# zones 


HLL 


HLLC 


Riemann 


# 1 


4000 


1 


1.05 


1.37 


#2 


4000 


1 


1.07 


1.34 


# 3 


4000 


1 


1.06 


1.30 


#4 


4000 


1 


1.04 


1.44 


R2-D 


400 2 


1 


1.03 


1.30 



Table 1. Normalized CPU time per integration step for the first 
four one-dimensional shock tubes and the 2-D Riemann problem 
(R2-D) considered in H4.7I The two rightmost columns give the 
average computing time for the HLLC and two-shock nonlinear 
Riemann solvers normalized to the HLL time (third column). All 
runs were produced using the first-order scheme with CFL = 0.8. 



fourth test problems is particularly time consuming, since a 
stronger discontinuity is involved. 

In this respect, a direct comparison between different 
Riemann solvers becomes problem-dependent and can be 
used only as an order-of-magnitude estimate. Conversely, 
we do not expect the HLLC/HLL efficiency ratio to change 
with increasing complexity of the flow patterns. For this rea- 
son, we believe that for problems involving rich and complex 
structures the trade-off between computational efficiency 
and quality of results is certainly worth the effort. 



6 CONCLUSIONS 

We have presente d, for the f i rst tim e, an extension of the 
HLLC scheme by iToro et alJ dl994h to relativistic gas dy- 
namics. The solver is robust, computationally efficient and 
complete, in that it considers the full wave structure in the 
solution to the Riemann problem. The solver retains the at- 
tractive feature of being positively conservative, typical of 
the HLL scheme family. 

The major improvement over the simple single-state 
HLL solver is the ability to resolve contact and tangential 
discontinuities. This property has been demonstrated by di- 
rect comparisons in several one- and two-dimensional test 
problems, where differences are strongly evident. The results 
indicate that the new HLLC solver attains sharper represen- 
tation of discontinuities, quantitatively similar to the exact 
but algebraically and computationally more intensive Rie- 
mann solver. 

The additional computational cost over the traditional 
HLL approach is less than 8% and we believe that the im- 
proved quality of results largely justifies the trade-off be- 
tween the two approximate Riemann solvers. 

Extension to relativistic magnetized flows will be con- 
sidered in a forthcoming paper. We notice, however, that 
the HLLC formalism presented in this work can be easily 
generalized to the case of vanishing normal component of 
magnetic field. When this degeneracy occurs (as in the prop- 
agation of jets with t oroidal magnetic field, see for example 
iLeismann et alj|2005l) . in fact, the solution to the Riemann 
problem is entirely analogous to the non-magnetized case, 
since only three wa ves are actually i nvolv ed. This extension 
will be presented in iMignone et alJfeOOSh . 



Finally, we mention that the relativistic HLLC scheme 
does not make any assumption on the equation of state, and 
efforts to incorporate different equations of state should be 
minimal. 
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APPENDIX A: PROOF OF Al < A* < Ar 

In what follows we prove some important results concern- 
ing the positivity of the our relativistic HLLC scheme. The 
proof is given below in IA.4I Propositions IA.1I through I A. 31 
demonstrate some preliminary results. 

We assume that the fastest and slowest signal velocities 
are computed according to the prescription given in 83.1.11 
and that Xr > and Al < 0, which is the case of applica- 
bility for the intermediate fluxes (1131 . Obviously, the initial 
left and right states are supposed to be physically admissi- 
ble, i.e. they belong to the set G defined in 83.1.21 

Proposition A.l. A R > 0, A L < 

Proof We will only prove Ar > 0, since the proof for Al 
is similar. For the sake of clarity, we omit the subscript R. 
Using the definition of A given after equation Ijl7^ one has 



A = (E + p) 



A 1 



(Al) 



where a = c 2 /^ 2 is always positive numbers. Equation lAH 
is always positive for A > Ao, where 



A = 



(A2) 



However, according to the choice given in 83.1.11 A must 
satify 



f(X) = (X-v x ) 2 



1 



r(l-A 2 ) 3*0. 



(A3) 



with a s = (t/(1 — cf). Equation 1A3I simply states that A 
must be greater than the root with the positive sign A+ 
(for the left state, Al is always less than the root with the 
negative sign A_). Direct substitution of Ao from <A2> into 
iA3l shows, after some algebra, that 



/(Ao) < Ki [-ci - 



2r) c j 



(A4) 



where the equality occurs in the limit of zero tangential ve- 
locities and K\ is always a positive quantity. Since 1 < T < 2 
and c 2 has the limiting value (r — 1) the expression in square 
bracket is always negative, which means that A ^ A+ > Ao. 
This implies that A = Ar is always positive with our choice 
of A = Ah. 

Since one can prove, in a similar way, that Al < 0, we 
have the important results that U^ u = Ar — Al > 0. 

Proposition A. 2. B R + A R > 0, B L - A L > 0. 

Proof Again we only give the proof for the right state, the 
other case being similar. The function A + B (the subscript 
R is omitted), with A and B defined after equation 11711 
increases linearly with A and is positive for A > Ao, where 



Ao = 



v x j 2 r(v x + 1) + c 2 s 



(A5) 



^T(v x + 1) - c; 

However, direct substitution of Ao in 1A3II shows, after ex- 
tensive manipulations, that 

/(Ao) = -K 2 



i + 2r)<^ + 7 2 r 2 (i- 



(A6) 



where Ki is always a positive quantity. It can be easily ver- 
ified that the function in square brackets is always posi- 
tive if c 2 s e [0, T - 1] and 1 < T < 2. Thus we must have 

Br/Ar > -1. 

Proposition A. 3. X L A R - Ar < 0, X R A L - B L > 0. 
Proof For the right state we have that 

X L A - B sC X-A - B sS ( -r^- x]A - B 



1 



(A7) 



where the last inequality follows from the fact that the two 
roots of equation JA3II satisfy 



A_ 



2v 



1 



— Aj 



and A+ sC A . 



(A8) 



Using the fact that A 2 (l + <x s ) ^ 2 Aw — v 2 + a s and that 
A > 0, the last expression in equation 1A7I can be shown to 
obey the following 

2v 



(—-A 

\1+(T S J 



.l + o- 
where 

9 = K 3 



-X)A-B^g, 



1 - r + 4 



(A9) 



(A10) 



with Ki being a positive quantity. The expression in square 
bracket in equation lAlOfl is always negative under the same 
assumptions used previously. Thus we have Al < Br/Ar 
and, similarly, one can prove that A_r > Bl/Al- 

Proposition A. 4. Al ^ A* ^ Xr. 

Proof We now show that the choice of eigenvalues given in 
83. 1.1 1 always guarantees Al ^ A* ^ Xr. 

The starting point is to note that the quadratic equation 
I18H can be more conveniently written as 

(AlA* - B L )(1 - AfjA*) = {ArX - B R )(1 - AlA*) , (All) 

which defines the intersection of two quadratic functions. 
The parabola on the left hand side vanishes in A* = 1/Xr > 
1 and A* = Bl/Al < I, whereas the parabola on the right 
hand side in A* = 1/Al < —1 and A* = Br/Ar > — 1. 
Moreover the two quadratics have the same concavity, since 
sign(ALAfi) = sign(AflAL) = — 1. Thus the intersection 
must necessarily satisfy 

B± 

A L , 

However, for any A G (— 1, 1) one has 
XA — B = (X — v x f{E + p) + p(l - A 2 ) > , 



(Br B L \ (Br 

— — , — — < A < max — — . 

\Ar A L ) \Ar' 



(A12) 



(A13) 

which, together with the results previously shown, implies 
that 



1 > A_r > max 



/ Br Bl 



\Ar' Al 



)• 



— 1 < Al < min 



(Br Bl 
\Ar' A l 



(A14) 
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and hence 

\l < A* < \ R . (A15) 



